CUBIT and Seismic Wave Propagation Based Upon the Spectral-Element Method: An Advanced Unstructured Mesher for Complex 3D Geological Media

Unstructured hexahedral mesh generation is a critical part of the modeling process in the Spectral-Element Method (SEM). We present some examples of seismic wave propagation in complex geological models, automatically meshed on a parallel machine based upon CUBIT (Sandia Laboratory, cubit.sandia.gov), an advanced 3D unstructured hexahedral mesh generator that offers new opportunities for seismologist to design, assess, and improve the quality of a mesh in terms of both geometrical and numerical accuracy. The main goal is to provide useful tools for understanding seismic phenomena due to surface topography and subsurface structures such as low wave-speed sedimentary basins. Our examples cover several typical geophysical problems: 1) “layer-cake” volumes with high-resolution topography and complex solidsolid interfaces (such as the Campi Flegrei Caldera Area in Italy), and 2) models with an embedded sedimentary basin (such as the Taipei basin in Taiwan or the Grenoble Valley in France). Emanuele Casarotti, Marco Stupazzini, Shiann Jong Lee, Dimitri Komatitsch, Antonio Piersanti, and Jeroen Tromp 580 Emanuele Casarotti et al. 1 The spectral-element method in seismology Wave propagation phenomena can nowadays be studied thanks to many powerful numerical techniques. We have seen rapid advances in computational seismology at global, regional, and local scales thanks to various numerical schemes, such as the finite-difference method (FDM) [15], the finite-element method (FEM) [1], the spectral-element method (SEM) [12], the integral boundary-element method (IBEM) [2], and the Arbitrary High-Order Discontinuous Galerkin method (ADER-DG) [8]. Spurred by the computational power made available by parallel computers, geoscientists and engineers can now accurately compute synthetic seismograms in realistic 3D Earth models. Among the methods previously listed, the SEM has convincingly demonstrated the ability to handle high-resolution simulations of seismic wave propagation in 3D domains. The purpose of this article is not to give a complete description of the SEM, but a basic introduction to the main properties of the SEM is needed in order to understand the various constraints imposed on the meshing process. The SEM is as a generalization of the FEM based on the use of high-order piecewise polynomial functions. A crucial aspect of the method is its capability to provide an arbitrary increase in accuracy by simply increasing the algebraic degree of these functions (the spectral degree n). From a practical perspective, this operation is completely transparent to the users, who limit themselves to choosing the spectral degree at runtime, leaving to the software the task of building up suitable quadrature points and the relevant degrees of freedom. Obviously, one can also use mesh refinement to improve the accuracy of the numerical solution, thus following the traditional finite-element approach. Referring to the literature for further details [12], we begin by briefly summarizing the key features of the SEM. We start from the general differential form of the equation of elastodynamics: ρ∂ t s = ∇ ·T+ F, (1) where s(x, t) denotes the displacement at position x and time t, ρ(x) is the density distribution, and F(x, t) the external body force. The stress tensor T is related to the strain tensor by the constitutive relation T = c : ∇s, where c denotes a fourth-order tensor. No particular assumptions are made regarding the structure of c, which describes the (an)elastic properties of the medium (the formulation is general and can incorporate full anisotropy [9] or non-linear constitutive relationships [6]). For seismological applications, a numerical technique needs to facilitate at least the following: (i) an attenuating medium, (ii) absorption of seismic energy on the fictitious boundaries of the domain in order to mimic a semiinfinite medium (the free-surface condition is a natural condition in the SEM), and finally (iii) seismic sources. In the SEM all these features are readily accommodated [12]. 5B.4 CUBIT and Seismic Wave Propagation 581 As in the FEM, the dynamic equilibrium problem for the medium can be written in a weak, or variational, form, and, through a suitable discretization procedure that depends on the numerical approach, can be cast as a system of ordinary differential equations with respect to time. Time-marching of this system of equations may be accomplished based upon an explicit second-order finite-difference scheme, which is conditionally stable and must satisfy the well known Courant condition [4]. The key features of the SEM discretization are as follows: 1. Like in the FEM, the model volume Ω is subdivided into a number of non-overlapping elements Ωe, e = 1, . . . , ne, such that Ω = ⋃ne e=1 Ωe. 2. The expansion of any function within the elements is accomplished based upon Lagrange polynomials of suitable degree n constructed from n + 1 interpolation nodes. 3. In each element, the interpolation nodes are chosen to be the GaussLobatto-Legendre (GLL) points, i.e., the n+1 roots of the first derivatives of the Legendre polynomial of degree n. On such nodes, the displacement, its spatial derivatives, and integrals encountered in the weak formulation are evaluated. 4. The spatial integration is performed based upon GLL quadrature (while most classical FEMs use Gauss quadrature). Thanks to this numerical strategy, exponential accuracy of the method is ensured and the computational effort minimized, because the resulting mass matrix is exactly diagonal. The latter feature does not occur in so-called hp FEMs, nor in SEMs based upon Chebychev polynomials. 2 Mesh design for Spectral-Element Methods The first critical ingredient required for a spectral-element simulation is a high-quality mesh appropriate for the 3D model of interest. This process generally requires discouraging expertise in meshing and preprocessing, and is subject to several tough constraints: 1) the number of grid points per shortest desired wavelength, 2) the numerical stability condition, 3) an acceptable distortion of the elements, 4) balancing of numerical cost and available computing resources. A poor-quality mesh can generate numerical problems that lead to an increase in the computational cost, poor (or lack of) convergence of the simulation, or inaccurate results. For example, a geological model often includes a layered volume, and a staircase sampling of the interfaces between the layers can produce fictitious diffractions. Therefore, a good mesh should honor at least the major geological discontinuities of the model. As noted in Section 1, the SEM is similar to a high-degree FEM, and in fact these methods share the first part of the discretization process. The present paper is focused on this first meshing step, and thus the results are relevant to both FEMs and SEMs. In the SEM, the subsequent step is the 582 Emanuele Casarotti et al. evaluation of the model on the GLL integration points [12]. Here we only note that if we use Lagrange polynomials and GLL quadrature, the mass matrix is exactly diagonal, resulting in a dramatic improvement in numerical accuracy and efficiency. With this choice, each element of the mesh contains (n + 1) GLL points. Unfortunately, this approach requires that the mesh elements are hexahedra [12]. It is worth mentioning that for 2D problems it is possible to develop a SEM with triangles keeping the diagonal mass matrix but with an higher numerical cost [10]. The fact that we are restricted to hexahedral elements complicates matters significantly. Whereas 3D unstructured tetrahedral meshes can be constructed relatively easily with commercial or non-commercial codes, the creation of 3D unstructured hexahedral meshes remains a challenging and unresolved problem. For complex models, as in the case of realistic geological volumes, generating an all-hexahedral mesh based upon the available meshing algorithms can require weeks or months, even for an expert user [21]. One of the features of the SEM that impacts the creation of the mesh is the polynomial degree n used to discretize the wave field. The following heuristic rule has emerged to select n for an unstructured mesh of a heterogeneous medium: if n < 4 the inaccuracies are similar to the standard FEM, while if n > 10 the accuracy improves but the numerical cost of the simulation becomes prohibitive. The choice of n is related to the grid spacing Δh: in order to resolve the wave field to a shortest period T0, the number of points per wavelength λ should be equal or greater than 5, leading to the constraint expressed in Eq. 2. If n = 4, then Δh is roughly equal to λ. We note that, for the same accuracy, a classical low-degree FEM requires a higher number of elements. Since the material properties are stored for each GLL point and can vary inside an element, we are able to interpolate the geological interfaces that our mesh is not able to follow. Nevertheless, this is an undesirable staircase sampling of the model, which introduces non-physical diffractions. Therefore, it is necessary that the mesh honors the major seismological contrasts. Furthermore, when a discontinuity is honored, the elements across the interface share some nodes, which will have the properties of the material below the discontinuity in one case, and the proprieties of the material above the discontinuity in the other case. Thus, the actual behavior of the seismic waves at the geological interface is perfectly mimicked in a way that cannot be achieved by an interpolation solely based upon Lagrange polynomials and the GLL quadrature. Another constraint on the design of the mesh is the stability condition imposed by the adoption of an explicit conditionally stable time scheme. For a given mesh, there is an upper limit on the time step above which the calculations are unstable. We define the Courant stability number of the time scheme as C = Δt(v/Δh)max, where Δt is the time step and (v/Δh)max denotes the maximum ratio of the compressional-wave speed and the grid spacing. The Courant stability condition ensures that calculations are stable if the Courant number is not higher than an upper limit Cmax [4], providing a condition that 5B.4 CUBIT and Seismic Wave Propagation 583 determines the time step Δt (Eq. 3). Again, an heuristic rule suggests that Cmax is roughly equal 0.3–0.4 for a deformed and heterogeneous mesh [12]. Like any technique based upon an unstructured deformed mesh, the SEM requires a smooth variation of the Jacobian and an acceptable distortion of the elements. Usually, to produce an acceptable accuracy the maximum equiangle skewness should not be greater than 0.5, although empirical tests show that s < 0.8 can sometimes be acceptable (Eq. 4) [12]. To sum up, spectral-element simulations require an unstructured allhexahedral conforming mesh subject to the following heuristic constraints: Δh = vmin T0 n+ 1 f(n) , (2) Δt < Cmax vmin vmax T0 n+ 1 f(n) , (3)